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ABSTRACT 

In this study, we have developed a computational method that allows numerical 
calculations of the time dependent compressible Navier-Stokes equations. The current 
results concern a study of flow past a semi-infinite flat plate. Flow develops from given 
inflow conditions upstream and passes over the flat plate to leave the computational 
domain without reflecting at the downstream boundary. Leading edge effects are 
included in this paper. In addition, specification of a heated region which gets con- 
vected with the flow is considered. The time history of this convection is obtained, 
and it exhibits a wave phenomena. 
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1. Introduction 

Computational methods applied to the Navier-Stokes equations have been a very important prob- 
lem in the field of aerodynamics. The study continues due to several unanswered aspects of specific 
flow situations, for example, compressible flow in open domains in the presence of a plate. 

The physical problem which initiated this study involved interaction of acoustic waves in the 
boundary layer. Such a phenomenon is important in the study of aerodynamic applications such as tur- 
boprop engine equipped aircrafts. They are known to have high fuel efficiency: unfortunately, they have 
high output of sound energy compared to aircrafts powered by conventional turbofan engines. Scien- 
tists working in the field wanted to understand if the noise from these turboprop engines has any 
favourable effects on the aircraft performance. In particular, do the sound waves interfere in a construc- 
tive sense, with the boundary layer formed during the flight? It is suspected that this phenomena is 
true. It is also suspected that the drag coefficient will change as the of the sound waves change. To 
answer this question, it is essential to have a fundamental study made available; unfortunately, there 
are difficulties. First these problems are hard to investigate in a mathematical sense. Secondly, it is 
difficult even to define acoustic sources which satisfy the field equations which are the set of perturbed 
compressible Navier-Stokes equations. Thus, instead of introducing acoustic sources in the flow, we 
propose to study generation of acoustic sources in some way, that may possibly interact with the boun- 
dary layer. Thus, we have addressed a simpler problem of a time dependent compressible fluid flow in 
which heating is imposed; however, no acoustic sources are introduced externally. 
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While these analyses appear to have been done in a sequence of papers by Rudy and Strikwerda 
[5,6], Bayliss, et al [2], and Abarbanel, et al [1], the results were actually constructed using a previously 
established boundary layer (using the code of Harris [4]) which defined the upstream or inflow boun- 
dary condition. In the current study, the flow, including the boundary layer development, is achieved 
by integrating the Navier-Stokes equations subject to appropriate upstream (inflow) boundary condi- 
tions. That is, the upstream values of density, temperature, and velocity components are prescribed in 
accordance with the physical conditions. At the downstream boundary of the computational domain, 
we eliminate any back flow using the nonreflective boundary condition that was developed by Rudy and 
Strikwerda [6]. If the flow is initialized by constant ambient values, then introducing a flat plate parallel 
to the flow yields the standard boundary layer solution as time increases. However, if one introduces 
different initial conditions in an isolated region, the flow tends to develop waves that pass through the 
computational domain as time increases. It is of interest to understand the nature of such waves. Our 
preliminary results indeed show existence of such waves. At the outset they can be thought of as 
acoustic waves. We plan to investigate the nature of such waves in a continuation of this study. The 
main purpose of this paper is to indicate the essential steps in constructing a code for this problem. 

2. Mathematical Formulation 

Let ft be the infinite region in the x-y plane. Let ft’ be the region that the plate occupies (see 
Fig. 1). If we denote the x and y components of the velocity by u(x,y,t) and v(x,y,t), the pressure by 
p(x,y,t), the density by p(x,y,t), the total energy by E(x,y,t) and the temperature by T(x,y,t), then we 
seek to solve the initial boundary value problem to determine 


U = [p,pu,pv,E] T in ft Ift' such that. 


au dF 3G 

3t + 3x + 3y 


Q 


( 2 . 1 ) 
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and where the stress tensor components and heat flux are given by 


2 3u 
-p +I tt ( 2- 


— ) 

dy 


T yy = 




3u . 
dx } 






= X 


y* 


Q« = 



Qy = 



The above equations are obtained when constant specific heats are assumed (calorically perfect gas). 
The boundary conditions are 


u = u„, v = 0, and p = p„ as x -> 


It - “ pc aT = ° as x °° ( 2 - 2 > 

and, 

u = v = 0 on d£2' (boundary of the plate). 

The initial conditions (at t = 0) are 

u = u„ , v = 0, p = p„ , T = T„ in Q . 

Equation (2.2) specifies the nonreflective condition in an asymptotic limit. This is obtained by 
considering the appropriate inflow Riemann invariants for the inviscid unidirectional flow. For no 
reflection, the inflow variants must be zero and we obtain (2.2). This procedure is due to Hedstr&m [8]. 
An improvement on this condition has been given by Rudy and Strikwerda [6]: 

It" “ pC '!t’ + 6(P_ P “ ) = °* (23) 

Here ’0’ is a parameter that needs to be chosen optimally and works effectively even if the dis- 
tance is finite. The boundary conditions (2.2) and (2.3) are in the nondimensional form; thus the 
governing equations (2.1) need to be cast in the same manner. This procedure is described in the 
Appendix A of this paper. The plan of this paper is as follows: 


i 
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Beginning with the nondimensional form of the equations and boundary conditions, we shall 
describe the numerical procedure in the next section. This will also include the treatment of the boun- 
dary conditions numerically. Finally, the results of computations will be presented. 

3. Numerical Procedure 

The numerical method, commonly known as MacCormack’s scheme, consists of an explicit finite 
difference method which is second order accurate in time and space. It is possible to increase the spa- 
tial accuracy of the method using an upgrade of the MacCormack’s scheme (Gottleib and Turkel [3]); 
however, the viscous terms that contain the mixed derivatives in both methods are second order accu- 
rate, and thus the fourth order accuracy can be achieved only in the inviscid region. Since the accuracy 
in the viscous region is the most important feature, we are confined only to the second order accurate 
scheme globally. 

For numerical computations, it is essential to truncate the infinite region Q to one of finite size 
which we will denote by Q . . An obvious choice is to truncate the region as in Figure 2. 

Let T { be the inflow boundary, r o be the outflow boundary, and let T sl and r s2 denote the upper 
and lower boundaries of the calculation domain respectively. If we assign the origin at the left hand bot- 
tom corner (as in Fig. 2), then the Navier-Stokes equations (2.1) are descretized according to the expli- 
cit two stage difference formula: 


U£j +1 = LLij- 


At 

Ax 


[Bij - EG ]- ^ [av. - Ofi ] - "H 


n 


(3.1) 


uy 1 = + m- [ay 1 - e^]- [sir - g&i]- 


(3.2) 


In the absence of body forces such as gravity, the nonhomogeneous terms will be zero. That is, 
H = 0. In these formulas. 


j-i 

Xj — (i - 1) Ax and yj = £Ay k 0< i< N, 0< j< M. 

k=0 

Thus, equations (3.1) and (3.2) require values of 


En+u » Si.M+1 > Eo.j and Gj >0 


which are obtained by extrapolating F and Q_ according to the relations (3.3). 
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En+ i,j - 2 F N j - En— ij 

^ (^Ym + AyM-l)Gi,M _ Ay M GiM-l 

iii,M+l - 77 r 

(Ay M -i) 

(3.3) 

Eoj - 2 Fjj - E 2 ,j 

0 (Ayo + Ay^Gi,! - Ay 0 Gi i2 

a '° " (*£> J 

It is important to note that the flux quantities F and Q. themselves contain derivatives with respect 
to x and y. Thus, the construction Ey and Crg need to be handled in a special way and they are given in 
Table 1. This is the same procedure given by Anderson, et al [9]. 


Predictor 

Ax 

Corrector 

fi+i.j - fij 
Ax 

Predictor 

fi+l.j~ fi— l,j 
2 Ax 

Corrector 

fi+l,j ~ fi-l,j 
2 Ax 


Table 1.x- derivatives where T denotes ’u,’ ’v,’ or ’T’ 
The y derivatives are handled globally as follows: 


3f 

Jy 


fj,j+l fj.j 
Ayj 

Ej ~ fj,j- 1 

Ayj-i 


.above the plate, 
.below the plate. 


The two differencing formulas in the y direction were motivated by the boundary conditions on 
the plate, which is discussed later. However, the accuracy is possibly decreased, and it requires further 
investigation. 

The above discussion completes the finite difference relations inside the computational domain. 
We also require the differencing relation on the physical and nonreflective boundary conditions men- 
tioned in section 2. 


4. Discrete Form of the Boundary Conditions 

Recall on the inflow boundary r ; , that p, u, and v were specified. In addition we specify T at the 
inflow. While this may be an over specification, this is one of the choices to make the initial value prob- 
lem well-posed. Details are seen in [6]. This translates into 



(4.1) 
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J (1 * jS M) 

In the set of equations (4.1), the quantities with subscripts °° denote the upstream values and are 
given in Appendix B. On the boundaries T sl and r j2 , the condition v = 0 is imposed to indicate the 
inviscid nature of the flow away from the plate. This, when written in discrete form, is 

vg =0 (1 < i< N) (j= 1 or M). (4.2) 

The outflow conditions are the same as that used by Rudy and Strikwerda [6]. In their work they 
show a set of four outflow boundary conditions against computational efficiency. In this work, the 
choice was the set which is mathematically most consistent but a bit slower in achieving convergence. 
Here we calculate the pressure from the discrete form of the nonreflective condition (2.3) and then cal- 
culate temperature using the equation of state for a gas which is calorically perfect (constant specific 
heats): 


Pt.j = P" 


1 i.j 
u i.j 
V U 


= T„ 
= u„ 
= 0 


p = (Y - l)pT 


(4.3) 


The discrete form of the nonreflective condition (2.3) is 


Pnj 1 - (j ~0At) + ® At P~ + PnjCnj ( u nj* ~ u N,j) j • (4.4) 

Using relation (4.3), the temperature is determined by 


T n+1 _ 

N,j - 


n n + 1 
Pnj 


(Y ~ 1)PnV ’ 

The speed of sound, C, is calculated according to (See Appendix A) 


(4.5) 


Cij = V Y(Y - l)Tij. 


(4.6) 


Finally u,v and p are prescribed by zeroth order extrapolation on the outflow boundary. 

Along the plate, both u and v are zero and the temperature of the plate is maintained at T„, . In 

addition it is necessary to impose — = 0 ( see [6] ). The discrete forms of these conditions are 

dy 


-7- 


Typ = T„ = T^jp-i 

Ui,jp = 0 = Uijp_i 
Vi,jP = 0 = Vyp.j 

Pi,jP+ 1 = Pi,jP 
Pi,jP-2 ~ PiJP- 1 


(4.7) 


where jP indicates the upper boundary of the plate and jP-1 is the lower one. 

The differencing used for y derivatives in the viscous terms is different than that suggested by 
Anderson, et al [9]. The reason we give for this change is as follows: the backward difference at jP or 
the forward difference at jP-1 yields 3u/3y = 0 which is not consistent with existing physical conditions 
on the plate (in fact, 3u/3y = 0 implies no velocity gradients at the plate which means in viscid flow!). 
To overcome this difficulty, we have used forward differencing in the y region above the plate and 
backward differencing in the region below. In addition to boundary conditions, initial conditions and 
reference values must be specified . 


5. Numerical Results and Discussions 

With the numerical scheme described in the last section, a code was developed. Calculation was 
started from a state of rest with initial conditions p^, u,*, v M , and T« . The solution was monitored for 
a change in solution of the order of 10" 6 in each dependent variable in the successive calculations. If 
is any of the dependent variables, then convergence was assumed when 


^ 10- 6 


To achieve this convergence, typically 60,000 steps were needed. It is also important to note that 
the Reynold’s number used in the calculation is 1.5 x 10 s for the configuration considered with 
L ref = 1 ft. 

At this point we remark that the time step At was chosen well below the stability limit that was 
prescribed by MacCormack [9]. A nondimensional time step of 10 -5 was chosen for the calculations 
for the meanflow and for the case in which there were no thermal disturbances. Moreover, to capture as 
many points as possible inside the boundary layer while keeping the upper boundary of the calculation 
domain as far away from the plate as possible, stretching was used. The stretching was obtained using 
the following formula. 


y = 



P 

P - i + y 


(5.1) 
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This stretching was above the plate and was reflected in the plate to get a symmetric grid. The 
plate was assumed to have a thickness of a mesh layer adjacent to it. The value of (3 used was 1.000001. 
Such a transformation is due to Roberts [10]. The total number of grid points in the x direction was 60 
and also 60 in the y direction. 

The first case here is to let the flow develop the boundary layer. Mach number for all cases con- 
sidered here was 0.4. The mesh thickness in the x direction was 0.3, and the computational domain has 
an aspect ratio of roughly .95. The tuning parameter at the outflow boundary condition 0 was chosen to 
be 1.1 and was found to be optimal by experimentation. Plate is assumed to have a mesh width in the y 
direction and is thin due to the stretching. The leading edge of the plate is located at the 9th grid point 
in the x direction. Prandtl number was 0.72. Uniform initial conditions were used in the entire compu- 
tational domain and the algorithm was initiated to reach a convergence of the order 10~ 6 in 60,000 time 
steps. The result was compared against the incompressible case Blasius’s solution and is given in Figure 
3. The comparison was made at the 25th x- grid point on the plate. We obtained 9 points in the boun- 
dary layer. Numerical values are given in Table 2. 

Since these computations contain the leading edge, and very little is known about this situation in 
the computational literature, it is interesting to see the flow field at the last step of the calculation. We 
chose to present the flow field in the entire computational domain. In Figure 4 we see the u velocity 
component. The valley of the graph is the boundary layer region. The flow is moving from left to right. 
The flanged level indictes the region where the value of u has reached 1, the uniform flow condition. 
At the leading edge, we see some disturbances as expected. At the trailing edge where we had imposed 
the nonreflective condition, we see some reflections, because the boundary condition is not completely 
nonreflective. This could be improved by a more accurate nonreflective condition, about which unfor- 
tunately very little is known in the literature. Alternatively, it could be improved by stretching the com- 
putational domain at farther distances which will increase the computational time substantially. Settling 
for the current results that we have, we took a closer look at the flow, particularly at the leading edge. 
We opted to look at the crossection through a contour plot, which is shown in Figure 5. Stagnation 
region is clearly seen from the figure, and the flow tends to spread at the leading edge. An even more 
interesting feature that we observed is that of the v (y - velocity ) component. This is shown in Figure 
6. As we see above and below the plate, the velocity seems to undergo a cycle in opposite directions at 
the leading edge. A reverse cycle seems to take place at the trailing edge. It suggests that a vorticity pat- 
tern is created at the leading edge and the trailing edge, where the nonreflective condition is imposed. 
Taking a closer look at the crossection through the contours, in Figure 7, we see periodic orbits 
corresponding to the disturbances at the leading edge. They also exist at the trailing edge, but exhibit 
inaccuracies in the nonreflective boundary condition. Such phenomena need further investigation, and 
we are unable to explain it fully in this paper. 
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To investigate the wave phenomena, we initialized a temperature which has a value twice that of 
the ambient value in a small region of 4 mesh size just above the leading edge. As the flow moves the 
temperature certainly gets convected along the flow and finally reaches a steady value. But in this transi- 
tion we suspect that the nature of convection should be of a wave phenomenon, in partcular an acoustic 
one. To understand this transition, we plotted the time history of the temperature, which we saved at 
every step of integration. However, it is known that the numerical scheme itself has its own nature of 
oscillations. To separate this effect, we plotted the time history of the temperature in the same region 
without initializing higher temperature. This is shown in Figure 8. We then computed the time history 
of the heated region, which is shown in Figure 9. It is clear that the time history in the heated region 
has a different frequency of oscillations than the one with no heat, exhibiting a true wave phenomenon. 
Finally we compared the case of the heated region and the case unheated region in the same graph 
(Figure 10). Clearly we see the oscillations that correspond to the numerical scheme and the ones that 
correspond to the physical situation. In particular, up to the value of time t = .1 , we see a high fre- 
quency oscillation. It should be noted that similar concepts are found in reference [7], where a simple 
driven cavity problem was considered for these purposes. 

It is curious to know what happens to the boundary layer in the presence of the heated region. 
After 60,000 time steps for each case (with and without heated region) we constructed the values of u 
- velocity component within the boundary layer as before. Although the numerical values do not differ 
by a substantial amount; there are differences. The results are given below. 
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TABLE 2 


Boundary Layer Solutions at i = 25 

Point 

Regular Flow 

Flow with Heat 

0 

0.00000 

0.00000 

1 

0.23222 

0.23221 

2 

0.45365 

0.45372 

3 

0.65929 

0.65933 

4 

0.81739 

0.81739 

5 

0.92020 

0.92023 

6 

0.97135 

0.97131 

7 

0.99170 

0.99168 

8 

0.99801 

0.99799 


* (0 represents the plate level) 

All our computations were performed on a SUN Microsystem (3/260) with a floating point 
accelerator. Double precision calculations took 65 CPU hours of computation in each case. 

From our computations, it is clear that two important things need to be investigated in the future. 
First, the leading edge singularity needs to be understood further by computational work with a theoret- 
ical backing. Second is of course, the nature of these waves; again to isolate the nature, theoretical 
ideas are needed. While there is a large amount of literature available to interpret these phenomena for 
incompressible flows, we are not aware of satisfactory theoretical work to explain these results for the 
compressible flow situation. We intend to continue this work in the future to answer these questions. 
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Appendix A 
Non dim en sionalization 


Let L be a characteristic length in the flow domain. Our coordinate distances x and y are nondi- 
mensionalized with respect to this length. We also define reference values for which all of the thermo- 
dynamic and fluid quantities are known. These are denoted by the subscript "ref," We note that these 
reference values and the conditions at negative infinity may very well be the same (but not necessarily). 
With this, we can write 


x = x/L y = y/L 


L/ u rc f 


u = u/u ref , V = v/v rcf 

jl = ft/ ft ref » P = p/Pref > P = P^Pref 

T = T/T re f , e = e/ e rC f , E = E/E re f . 


Some of the above quantities are not independent. The relationships are as follows: 

C ref = tl rc f 


2 

Pref “ Pref ^ref 


E re f — Pref 6 re f — Pref ^ref • 

The equation of state, when subjected to the above nondimensionalization, gets converted as fol- 
lows: 

p = p R T 


P = 


Pref R T ref 


pT. 


Pref 
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Dropping the bar notation for convenience. 


P = P PT 


where, 


PrefRT ref 

Pref 


RT 


ref 


U r ef 


Now, 

ui 

T ref = — and u ref = u„ 
c v 

(Note that we have equated the reference plane and infinity) which yields 

P = (Y - 1) 


and 


P = (Y " 1)PT. 


The Navier-Stokes equation takes the following nondimensional form: 


3U 

at 


+ 




( 2 . 1 ) 


where 


U = 


P 

pu 

pv 
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pu 

pu 2 — x xx 
puv - t xy 
_ 3 T 

E U - Pa-^- - UT XX - VT xy 


G = 


pv 

pUV-Ty, 

PV 2 - Tyy 

_ 3T 

Ev-p a — - UTy X - VX xy 


and where the stress tensor components and energy are given by 
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E 


e + 


u 2 + v 2 




X 


yy 


-P + 


-P + 


2 |I /pdii _ 3v 

3 Re 3x 3y 

2 |x 3v _ 3u 

3 Re dy dx 


‘'xy 


= = x 

Re 3y + 3x 


yx 


and 


Ha 


c p T rcf 
Pr Re u£ f 


n Pref u refL ^ CpHref 

Re = , rT— : . 


Href 


^ref 


In these calculations, we have used \ia = p,. 
To calculate the speed of sound, ’C,* 


C = 


11/2 


TP 


[Xll ^P T ] 1/2 since p = (y - l)pT 

P 


= V Y(y- 1)T. 
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